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f"*) ' Abstract. 

^SJ , Feedback from photoionisation may dominate on parsec scales in massive star-forming re- 

gions. Such feedback may inhibit or enhance the star formation efficiency and sustain or even 
O drive turbulence in the parent molecular cloud. Photoionisation feedback may also provide a 
, mechanism for the rapid expulsion of gas from young clusters' potentials, often invoked as the 
main cause of 'infant mortality'. There is currently no agreement, however, with regards to the 
efficiency of this process and how environment may affect the direction (positive or negative) 
in which it proceeds. The study of the photoionisation process as part of hydrodynamical simu- 
I I ■ lations is key to understanding these issues, however, due to the computational demand of the 
\ problem, crude approximations for the radiation transfer are often employed. 

OWe will briefly review some of the most commonly used approximations and discuss their ma- 
jor drawbacks. We will then present the results of detailed tests carried out using the detailed 
photoionisation code MOCASSIN and the SPH-fionisation code iVINE code, aimed at under- 
CIh' standing the error introduced by the simplified photoionisation algorithms. This is particularly 
1^ ' relevant as a number of new codes have recently been developed along those lines. 
^ , We will finally propose a new approach that should allow to efficiently and self-consistently 



treat the photoionisation problem for complex radiation and density fields. 
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Tj" 1. Introduction 

Ionising radiation from OB stars influences the surrounding interstellar medium (ISM) 
ly-j on parsec scales. As the gas surrounding a high mass star is heated, it expands forming 
an HII region. The consequence of this expansion is twofold, on the one hand gas is 
I removed from the centre of the potential, preventing further gravitational collapse and 
perhaps even disrupting the parent molecular cloud. On the other hand gas is swept up 
T-H and compressed beyond the ionisation front producing high density regions that may be 
^ ' susceptible to gravitation collapse (i.e. the "collect and collapse" model, Elmegreen et 
• 'H al. 1995). Furthermore, pre-existing, marginally gravitationally stable clouds may also 
K% be driven to collapse by the advancing ionisation front (i.e. "radiation-driven implosion" , 
^ . Bertoldi 1989). Finally, ionisation radiation has also been suggested as a driver for small 
■ " " ' scale turbulence in a cloud (Gritschneder et al. 2009b). Observations (e.g. Deharveng, 
these proceedings) and theory (e.g. Dale et al. 2005, 2007, Gritschneder et al. 2009b) 
often present examples for positive and negative feedback, however, the net effect on the 
global star formation efficiency is still under debate. 

From a theoretical point of view, different groups have performed a number of numeri- 
cal experiments demonstrating that the efficacy and direction of photoionisation feedback 
are very sensitive to the specific initial conditions, in particular, to the location of the 
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ionising source(s) and to whether the cloud is initially bound or unbound. This suggests 
that a parameter space study may be necessary to assess what environmental variables 
may affect the direction in which feedback proceeds. Several authors in these proceed- 
ings discuss the results of recent ionisation feedback simulations (see oral contributions 
by Arthur, Bisbas, Gritschneder and Walch, and poster contributions by Ghoudhury, 
Cornwall, A^Eiao, Motoyama, Rodon and Tremblin). 

As the field matures and the codes become more sophisticated it becomes important 
to assess the accuracy and limitations of the methods currently employed. The computa- 
tional demand of treating the radiation transfer (RT) and photoionisation (PI) problem 
within a large scale hydrodynamical simulation has led to the development of approxi- 
mate algorithms that drastically simplify the physics of RT and PI. In this review we will 
describe some of the most common approximations employed by current RT+PI imple- 
mentations, highlighting some potentially important shortcomings. We will then present 
the result of our ongoing efforts to test current implementations against the 3D Monte 
Carlo code MOCASSIN (Ercolano et al. 2003, 2005, 2008) which includes all the necessary 
micro physics and solves the ionisation, thermal and statistical equilibrium in detail. 

2. Some Common Approximations 

The importance of studying the photoionisation process as part of hydrodynamical star 

formation simulations has long been recognised. Until very recently, however, due to the 
complexity and the computational demand of the problem, the evolution of ionised gas 
regions had only been studied in rather idealised systems (e.g. Yorke et al. 1989; Garcia- 
Scgura & Franco 1996), with simulations often lacking resolution and dimensions. The 
situation in the latest years has been rapidly improving, however, with more sophisticated 
implementations of ionised radiation both in grid-based codes (e.g. Mellcma et al. 2006; 
Peters et al. 2010) and Smoothed Particle Hydrodynamical (SPH) codes (e.g. Kcssel- 
Deynet & Burkert 2000; Miao et al. 2006; Dale et al. 2007; Gritschneder et al. 2009; 
Bisbas et al. 2009). Klessen et al. (2009) and Mac Low et al. (2007) present recent 
reviews of the numerical methods employed. 

While the new codes can achieve higher resolutions and can treat more realistic geome- 
tries, the treatment of RT and PI is still rather crude in most cases. Even in the current 
era of parallel computing, an exact solution of the radiative transfer (RT) and photoioni- 
sation (PI) problem in three dimensions within SPH calculations is still prohibitive. Some 
common approximations include the following: 

(a) Monochromatic radiation field: In order to avoid the burden of frequency resolved 
RT calculations, monochromatic calculations are often carried out, where all the ionising 
flux is assumed to be at 13.6 eV (i.e. the H ionisation potential). This approximation 
is often implicit in the choice of a single value for the gas opacity, and it is of course 
implicit to Stromgren-type calculations. Implicit or explicit monochromatic fields have 
the serious drawback that the ionisation and temperature structure of the gas cannot be 
calculated. 

(6) Ionisation and thermal balances: Its equations are not solved or simple heat- 
ing/cooling functions are employed or the temperature is a simple function of an approx- 
imate ionisation fraction. When monochromatic fields are employed it is not possible to 
calculate the necessary terms to solve the balance equations and idealised temperature 
distributions must be used. 

(c) On-the-spot (GTS) approximation (no diffuse field): The OTS approximation is 
described in detail by Osterbrock & Ferland (2006, page 24). In the OTS approximation 
the diffuse component of the radiation field is ignored under the assumption that any 



lonisation Feedback 



121 



ionising photon emitted by the gas wiU be reabsorbed elsewhere, close to where it was 
emitted, hence not contributing to the net ionisation of the nebula. This is not a bad 
approximation in the case of reasonably dense homogeneous or smoothly varying density 
fields, but it is certain to fail in the highly inhomogeneous star-forming gas, where the 
ionisation and temperature structure of regions that lie behind high density clumps and 
filaments is often dominated by the diffuse field. 

{d) Steady-state calculations (instantaneous ionisation): The ionisation structure and 
the gas temperature of a photoionised region is often obtained by simultaneously solving 
the steady state thermal balance and ionisation cqiiilibrium equations. This approxima- 
tion is valid when the atomic physics timescales are shorter than the dynamical timescales 
and the rate of change of the ionising field. In this case, the photoionisation problem is 
completely decoupled from the dynamics and it can be solved for a given gas density 
distribution obtained as a snapshot at a given time in the evolution of a cloud. This is a 
fair assumption for the purpose to study of ionisation feedback on large scales, as most 
of the gas will be in equilibrium. Non-equilibrium effects, however, should still be kept 
in mind when interpreting the spectra of regions close to the ionisation front or where 
shocks are present. 

3. How good are the approximations? 

In cases where the steady-state calculations are relevant, it is possible to test the effects 
of approximations a-c from the above list by comparing the temperature distributions 
obtained by the hydro-|-ionisation codes against those obtained by a specialised pho- 
toionisation code, like the MOCASSIN code, for density snapshots at several times in the 
hydrodynamics simulations. 

MOCASSIN is a fully three-dimensional photoionisation and dust radiative transfer code 
that employs a Monte Carlo approach to the fully frequency resolved transfer of radia- 
tion. The code includes all the microphysical processes that influence the gas ionisation 
balance and thermal balance as well as those that couple the gas and dust phases. In 
the case of an HII region ionised by an OB star the dominant heating process for typical 
gas abundances is H photoionisation, balanced by cooling via collisionally excited line 
emission (dominant), recombination line emission and free-bound and free-free emission. 
The atomic database included in MOCASSIN includes opacity data from Verner et al. 
(1996), energy levels, collision strengths and transition probabilities from Version 5.2 of 
the CHIANTI database (Landi et al. 2006, and references therein) and the improved 
hydrogen and helium free-bound continuous emission data of Ercolano & Storey (2006). 

Dale et al. (2007, DEC07) performed detailed comparisons against mocassin's solution 
for the temperature structure of a complex density field ionised by a newly born massive 
star located at the convergence of high density accretion streams. They found that the 
two codes were in fair agreement on the ionised mass fractions in high density regions, 
while low density regions proved problematic for the DEC07 algorithm. The temperature 
structure, however, was poorly reproduced by the DEC07 algorithm, highlighting the 
need for more realistic prescriptions. For more details see DEC07. 

More recently we have used the MOCASSIN code to calculate the temperature and 
ionisation structure of the turbulent ISM density fields presented by Gritschneder et al. 
(2009b, hereafter: G09b). The SPH particle fields were obtained with the iVine code 
(Gritschneder et al. 2009a) and mapped onto a regular 128^ Cartesian grid. In order 
to compare with iViNE, which calculates the RT along parallel rays, the stellar field 
in MOCASSIN was forced to be plane parallel, while the following RT was performed in 
three dimensions thus allowing for an adequate representation of the diffuse field. The 



122 



Barbara Ercolano & Matthias Gritschneder 




Figure 1. Surface density of electrons projected in the z-direction for the G09b turbulent ISM 
simulation at t = 0.5 Myr. Left: iVine; Middle: MOCASSIN H-only; Right: MOCASSIN nebular 
abundances. 

incoming stellar field was set to the value used by G09b {Q'jj — 5x10^ ionising photons 
per second) and a blackbody spectrum of 40kK was assumed. We run H-only simulations 
(referred to as "H-only") and simulations with typical HH region abundances (referred 
to as "Metals"). The elemental abundance are as follows, given as number density with 
respect to Hydrogen: He/H = 0.1, C/H = 2.2e-4, N/H = 4.0e-5, 0/H = 3.3e-4, Ne/H = 
5.0e-5, S/H = 9.0e-6. 

The resulting MOCASSIN temperature and ionisation structure grids were compared to 
those obtained by iViNE in order to address the following questions: 

(a) Are the global ionisation fractions accurate? 

(b) How accurate is the gas temperature distribution? 

(c) What is the effect of the diffuse field? 

(d) How can the algorithm be improved? 

3.1. Global Properties 

Figure 1 shows the surface density of electrons projected in the z-direction for the G09b 
turbulent ISM simulation at t = 0.5 Myr. The figure shows that no significant differences 
are noticeable in the integrated ionisation structure, implying that the global ionisation 
structure is correctly determined by iVine. This is also confirmed by the comparison of 
the total ionised mass fractions: at t = 0.5 Myr, iVine obtains a total ionised mass of 
13.9%, while MOCASSIN "H-only" and "Metals" obtain 15.6% and 14.0%, respectively. 
The agreement at other time snapshots is equally good (e.g. at t = 250kyr iVine obtains 
9.1% and MOCASSIN "Metals" 9.5%). 

It may at first appear curious that the agreement should be better between iVlNE 
and MOCASSIN "Metals" , rather than MOCASSIN "H-only" , given that only H-ionisation 
is considered in iVine. This is however simply explained by the fact that iVine adopts a 
"ionised gas temperature" [Thot] of lOkK, which is close to a typical HH region temper- 
ature, with typical gas abundances. The removal of metals in the "H-only" simulations 
causes the temperature to rise to values close to 17kK, due to the fact that cooling be- 
comes much less efficient without coUisionally excited lines of oxygen, carbon etc. The 
hotter temperatures in the "H-only" models directly translate to slower recombinations, 
as the recombination coefficient is proportional to the inverse square root of the tem- 




Figure 2. Density and temperature maps for the z = 25 slice of the G09 turbulent ISM simula- 
tion at t = 0.5 Myr. Top left: Gas density map; Top right: electron temperature, Te as calculated 
by iVine; Bottom left: electron temperature, as calculated by MOCASSIN with H-only; Bottom 
right: electron temperature, Te as calculated by MOCASSIN with nebular abundances. 



perature. As a result of slower recombination the "H-only" grids have a slightly larger 
ionisation degree. 

3.2. Ionisation and temperature structure 

Accurate gas temperatures are of prime importance as this is how feedback from ionising 
radiation impacts on the hydrodynamics of the system. In Figure 2 we compare the 
electron temperatures Tg calculated by iVine and MOCASSIN ("H-only" and "Metals") 
in a z-slice of the t — 0.5 Myr grid. The top-right panel shows the number density [cm~^] 
map for the selected slice. The large shadow regions behind the high density clumps 
are immediately evident from both figures. These shadows are largely reduced in the 
MOCASSIN calculations as a result of diffuse field ionisation. The diffuse field is softer than 
the stellar field and therefore temperatures in the shadow regions are lower. The higher 
temperatures in the shadow regions of the MOCASSIN "Metals" model are a consequence 
of the Helium Lyman radiation and the heavy elements free-bound contribution to the 
diffuse field. The rise in gas temperature shown in the MOCASSIN results at larger distances 
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Figure 3. Density slice at 250 kyr for the OTS iVine (left) and the diffuse field iVine (right). 
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Figure 4. Turbulence spectrum obtained for the standard OTS I Vine (solid lines), the control 
run with no ionising radiation (dotted line) and the diffuse field I Vine (dashed lines). 

from the star is not surprising and a simple consequence of radiation hardening and the 
recombining of some of the dominant cooling ions. 

4. Towards more realistic algorithms 

As iViNE solves the transfer along plane parallel rays, it has currently no means of 
bringing ionisation (and hence heating) to regions that lie behind high density clumps. 
This creates a large temperature (pressure) gradient between neighbouring direct and 
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diffuse- field dominated regions, which may have important imphcations for the dynamics, 
particularly with respect to turbulence calculations. The same problem is faced by all 
codes that employ the OTS approximation and thus ignore the diffuse field contributions. 

In order to investigate whether the error introduced by OTS approximation actually 
bears any consequence on the dynamical evolution of the system and on the turbulence 
spectrum, we propose here a simple zeroth order strategy to include the effects of the 
diffuse field in iVine and which can be readily extended to other codes. It consists of 
the following steps: (i) identify the diffuse field dominated regions (shadow); (ii) study 
the realistic temperature distribution in the shadow region using fully frequency resolved 
three-dimensional photoionisation calculations performed with MOCASSIN and parame- 
terise the gas temperature in the shadow regions as a function of (e.g.) gas density; (iii) 
implement the temperature parameterisation in iVlNE and update the gas temperatures 
in the shadow regions at every dynamical time step accordingly. 

We note that this approach allows for environmental variables, such as the hardness 
of the stellar field and the metallicity of the gas to be accounted for in the SPH calcula- 
tion, since their effect on the temperature distribution is folded in the parameterisation 
obtained with mocassin. 

Figure 3 shows a slice of the density structure snapshots at 250k year for a standard 
IVINE (left panel) compared to a first attempt at a diffuse field implementation in iVINE 
(right panel). These preliminary results indicate that the diffuse field affects the evolution 
of structures, promoting the detachment of clumps from the pillars, presumably by ex- 
cavating them from behind.. Slightly higher densities are achieved in some of the clumps 
in the diffuse field simulation at this age, suggesting that the formation of stars may 
be thus accelerated. The turbulence spectrum obtained in simiilations with or without 
diffuse fields arc rather similar, as shown in Figure 4, where the specific kinetic energy is 
plotted as a function of wave number in the case of the control run with no ionisation at 
aU (dotted lines), OTS iVine (solid lines) and diffuse field iVlNE (dashed fines). There 
is however tentative evidence for less efficient driving at the smaller scales, due to the 
fact that the large temperature gradients created by the OTS at the shadow regions are 
removed when the diffuse field is considered. 

We stress that the results presented here are to be considered only a first exploratory 
step to establish whether diffuse field effects are likely to play a role in the dynami- 
cal evolution of a turbulent medium. More detailed comparisons will be presented in a 
forthcoming article (Ercolano & Gritschneder 2010, in prep) 

5. Conclusions 

We have presented a review of the current implementations of photoionisation al- 
gorithms in star formation hydrodynamical simulation, highlighting some of the most 
common approximation that are employed in order to simplify the radiative transfer and 
photoionisation problems. 

We discuss the robustness of the temperature fields obtained by such methods in light 
of recent tests against detailed 3D photoionisation calculations for complex density dis- 
tributions typical of star forming regions. We conclude that while the global ionised mass 
fractions obtained by the simplified methods are roughly in agreement, the temperature 
fields are poorly represented. In particular, the assumption of the OTS approximation 
may lead to unrealistic shadow regions and extreme temperature gradients that affect 
the dynamical evolution of the system and, to a lesser extent, its turbulence spectrum. 

We propose a simple strategy to provide a more realistic description of the tempera- 
ture distribution based on parameterisations obtained with a dedicated photoionisation 
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code, MOCASSIN, which includes frequency resolved 3D radiative transfer and all the mi- 
crophysical process needed for an accurate calculation of the temperature distribution 
of ionised regions. This computationally inexpensive method allows to include the ther- 
mal effects of diffuse field, as well as accounting for environmental variables, such as gas 
metallicity and stellar spectra hardness. 
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